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ABSTRACT 

We present the 8th Pull Pocal Plane simulation set (PPP8), deployed in support of the Planck 2015 results. PPP8 consists of 10 fiducial mission 
realizations reduced to 18 144 maps, together with the most massive suite of Monte Carlo realizations of instrument noise and CMB ever generated, 
comprising 10"^ mission realizations reduced to about 10^ maps. The resulting maps incorporate the dominant instrumental, scanning, and data 
analysis effects; remaining subdominant effects will be included in future updates. Generated at a cost of some 25 million CPU-hours spread across 
multiple high-performance-computing (HPC) platforms, PPP8 is used for the validation and verification of analysis algorithms, as well as their 
implementations, and for removing biases from and quantifying uncertainties in the results of analyses of the real data. 

Key words, cosmology: cosmic background radiation - cosmology: observations - methods: data analysis -methods: high performance computing 


1. Introduction 

Planck^ is the third satellite to study the cosmic microwave back¬ 
ground (CMB) radiation. 

Launched in May 2009, Planck started its science observa¬ 
tions from the L 2 Lagrange point in August 2009 and completed 


^ Corresponding author: J. Borrill jdborrill@lbl. gov 
^ Planck (http://www.esa.int/Planck) is a project of the 
European Space Agency (ESA) with instruments provided by two sci¬ 
entific consortia funded by ESA member states and led by Principal 
Investigators from Rrance and Italy, telescope reflectors provided 
through a collaboration between ESA and a scientific consortium led 
and funded by Denmark, and additional contributions from NASA 
(USA). 


an all-sky survey (Planck Collaboration I 2014) approximately 
every six months until it was decommissioned in October 2013. 
Planck carried two instruments. The High Erequency Instrument 
(HEI; Lamarre et al. 2010; Planck HEI Core Team 2011), com¬ 
prising 52 detectors^ at six frequencies (100, 143, 217, 353, 545, 
and 857 GHz), completed its observations when its cryogens 
were exhausted in January 2012 after almost five surveys, while 
the Low Erequency Instrument (LEI; Bersanelli et al. 2010; 
Mennella et al. 2011), comprising 22 detectors at three frequen¬ 
cies (30, 44, and 70 GHz) continued operating throughout the 
satellite lifetime, completing more than eight surveys. 


^ The data from two detectors proved to be unusable, so only 50 de¬ 
tectors are included in the analysis. 
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The second release of Planck data (hereafter PR2-2015) is 
based on five HFI and eight LFI surveys (LFI Survey 9 having 
been reserved for atypical scans designed to assist in the under¬ 
standing of systematic effects), and for the first time includes 
results in polarization as well as intensity. PR2-2015 is accom¬ 
panied by a suite of papers, of which this is one, together with 
an online explanatory supplement, an ESA legacy data archive 
(http://pla.esac.esa.int/pla/) and NASA partial mir¬ 
ror (http://irsa.ipac.caltech.edu/Missions/planck. 
html), as well as resources to access and manipulate the 
full simulation suite described here (http://crd.lbl.gov/ 
cmb-data). 

Simulations play a number of important roles in the analysis 
of these Planck data, including: 

1. validating and verifying tools used to measure instrument 
characteristics by simulating data with known instrument 
characteristics, applying the tools used on the real data to 
measure these, and verifying the accuracy of their recovery; 

2. quantifying systematic effect residuals by simulating data 
with some particular systematic effect included, applying the 
treatment used on the real data to ameliorate that effect, and 
measuring the residuals; 

3. validating and verifying data analysis algorithms and their 
implementations by simulating data with known science in¬ 
puts (cosmology and foreground sky) and detector charac¬ 
teristics (beam, bandpass, and noise spectrum), applying the 
analyses used on the real data to extract this science, and 
verifying the accuracy of its recovery; 

4. debiasing and quantifying uncertainties in the analysis of the 
real data by generating massive sets of Monte Carlo (MC) 
realizations of both the noise and the CMB and passing them 
separately through the analyses used on the real data to quan¬ 
tify biases and uncertainties. 

The first two items are instrument-specific, and distinct 
pipelines have been developed and employed by the LFI and 
HFI Data Processing Centres (DPCs); details are provided in 
the DPC processsing papers (Planck Collaboration II 2015; 
Planck Collaboration III 2015; Planck Collaboration IV 2015; 
Planck Collaboration V 2015; Planck Collaboration VI 2015; 
Planck Collaboration VII 2015; Planck Collaboration VIII 
2015). However, the last two items require consistent simula¬ 
tions of LFI and HFI in tandem; such simulations are beyond the 
scope of either single-instrument pipeline. Furthermore, gener¬ 
ating the Monte Carlo simulations is the most computationally 
intensive part of the Planck data analysis and requires computa¬ 
tional capacity and capability far beyond those available at the 
DPCs. As a result, a massively parallel cross-instrument suite of 
codes has been developed to provide the collaboration with self- 
consistent simulations of all detectors at all frequencies - ref- 
ered to as the Full Focal Plane (FFP) - primarily using the high 
performance computing (HPC) resources at the National Energy 
Research Scientific Computing Center^ (NERSC) in the USA 
and at CSC-IT Center for Science"^ (CSC) in Finland. Since its 
first deployment in 2006, this suite of codes has been used to 
generate a sequence of FFP simulation-sets of increasing verac¬ 
ity, complexity, and volume. The second Planck data release is 
supported by FFP8. 

This paper is laid out as follows. Section 2 outlines the over¬ 
all specification of FFP8, Sect. 3 details the inputs to FFP8, 
Sect. 4 describes the three major pipelines used to generate 

^ http://www.nersc.gov 
http://www.csc.fi 


FFP8, Sect. 5 covers the key HPC aspects of FFP8, and Sect. 6 
summarizes the results. 

2. Specifications 

The FFP8 simulations include a set of fiducial mission realiza¬ 
tions together with separate sets of Monte Carlo realizations 
of the CMB and the instrument noise. They contain the domi¬ 
nant instrumental (detector beam, bandpass, and correlated noise 
properties), scanning (pointing and fiags), and analysis (map¬ 
making algorithm and implmentation) effects. 

In addition to the baseline maps made from the data from all 
detectors at a given frequency for the entire mission, there are a 
number of data cuts that are mapped both for systematics tests 
and to support cross-spectral analyses. These include: 

- detector subsets (“detsets”), comprising the individual unpo¬ 
larized detectors and the polarized detector quadruplets cor¬ 
responding to each leading/trailing horn pair;^ 

- mission subsets, comprising the surveys, years, and half¬ 
missions, with exact boundary definitions given in Planck 
Collaboration II (2015) and Planck Collaboration VII (2015) 
for LFI and HFI, respectively; and 

- half-ring subsets, comprising the data from either the first or 
the second half of each pointing-period ring. 

The various combinations of these data cuts then define 1 134 
maps, as enumerated in the top section of Table 1. The differ¬ 
ent types of map are then named according to their included de¬ 
tectors (channel or detset), interval (mission, half-mission, year 
or survey), and ring-content (full or half-ring); for example the 
baseline maps are described as channel/mission/full, etc. 

The fiducial realizations include instrument noise 
(Sect. 3.1.2), astrophysical foregrounds (Sect. 3.2), and the 
lensed scalar, tensor, and non-Gaussian CMB components 
(Sect. 3.3), and are primarily designed to support the valida¬ 
tion and verification of analysis codes. To test our ability to 
detect tensor modes and non-Gaussianity, we generate five 
CMB realizations with various cosmologically interesting - 
but undeclared - values of the tensor-to-scalar ratio r and 
non-Gaussianity parameter fyi (Table 2). To investigate the 
impact of differences in the bandpasses of the detectors at any 
given frequency, the foreground sky is simulated using both the 
individual detector bandpasses and a common average band¬ 
pass, to include and exclude the effects of bandpass mismatch. 
To check that the PR2-2015 results are not sensitive to the 
exact cosmological parameters used in FFP8 we subsequently 
generated FFP8.1, exactly matching the PR2-2015 cosmology. 

Since mapmaking is a linear operation, the easiest way to 
generate all of these different realizations is to build the full set 
of 1 134 maps of each of six components: 

- the lensed scalar CMB (cmb_scl); 

- the tensor CMB (cmb_ten); 

- the non-Gaussian complement CMB (cmbmgc); 

- the forgreounds including bandpass mismatch (fg_bpm); 

- the foregrounds excluding bandpass mismatch (fg jnobpm); 

- the noise (noise). 

We then sum these, weighting the tensor and non-Gaussian com¬ 
plement maps with ^fr and /nl, respectively, and including one 
of the two foreground maps, to produce 10 total maps of each 
type. The complete fiducial data set then comprises 18 144 maps. 

^ HFI sometimes refers to full channels as detsetO; here detset only 
refers to subsets of detectors. 
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Table 1. The numbers of fiducial, MC noise and MC CMB maps at each frequency by detector subset, data interval, and data cut. 


FFP8 and FFP8.1 Fiducial Maps 

Detectors 

Interval 

Cut 

Frequency 

Total 

30 

44 

70 

100 

143 

217 

353 

545 

857 






Mission 

Full 

1 

1 

1 

1 

1 

1 

1 

1 

1 

9 



Half-Ring 

2 

2 

2 

2 

2 

2 

2 

2 

2 

18 


Half-Mission 

Full 

2 

2 

2 

2 

2 

2 

2 

2 

2 

18 

Channel 


Half-Ring 

4 

4 

4 

4 

4 

4 

4 

4 

4 

36 

Year 

Full 

4 

4 

4 

2 

2 

2 

2 

2 

2 

24 



Half-Ring 

8 

8 

8 

4 

4 

4 

4 

4 

4 

48 


Survey 

Full 

8 

8 

8 

4 

4 

4 

4 

4 

4 

48 


Half-Ring 

16 

16 

16 

8 

8 

8 

8 

8 

8 

96 


Mission 

Full 



3 

2 

5 

6 

6 

3 

4 

29 



Half-Ring 



6 

4 

10 

12 

12 

6 

8 

58 


Half-Mission 

Full 



6 

4 

10 

12 

12 

6 

8 

58 

DetSet 


Half-Ring 



12 

8 

20 

24 

24 

12 

16 

116 

Year 

Full 



12 

4 

10 

12 

12 

6 

8 

64 


Half-Ring 



24 

8 

20 

24 

24 

12 

16 

128 


Survey 

Full 



24 

8 

20 

24 

24 

12 

16 

128 


Half-Ring 



48 

16 

40 

48 

48 

24 

32 

256 

Total 

45 

45 

180 

81 

162 

189 

189 

108 

135 

1134 


FFP8 Noise MC Maps 

Detectors 

Interval 

Cut 

Frequency 

Total 

30 

44 

70 

100 

143 

217 

353 

545 

857 

Channel 

Mission 

Full 

Half-Ring 

10000 

20000 

10000 

20000 

10000 

20000 

10000 

200 

10000 

200 

10000 

200 

10000 

200 

10000 

200 

10000 

200 

90000 

61200 

Half-Mission 

Full 

Half-Ring 

2000 

4000 

2000 

4000 

2000 

4000 

2000 

400 

2000 

400 

2000 

400 

2000 

400 

2000 

400 

2000 

400 

18 000 

14400 

Year 

Full 

Half-Ring 

4000 

8 000 

4000 

8 000 

4000 

8 000 

2000 

400 

2000 

400 

2000 

400 

2000 

400 

2000 

400 

2000 

400 

24000 

26400 

Survey 

Full 

Half-Ring 

8 000 

16000 

8 000 

16000 

8 000 

16000 

4000 

800 

4000 

800 

4000 

800 

4000 

800 

4000 

800 

4000 

800 

48 000 

52 800 

DetSet 

Mission 

Full 

Half-Ring 



3 000 

6000 

2000 

400 

5 000 

1000 

6000 

1200 

6000 

1200 

3 000 

600 

13 000 

800 

38 000 

11200 

Half-Mission 

Full 

Half-Ring 



1800 

3 600 

4000 

800 

10000 

2000 

12000 

2400 

12000 

2400 

6000 

1200 

8 000 

1600 

53 800 

14000 

Year 

Full 

Half-Ring 



3 600 

7 200 

4000 

800 

10000 

2000 

12000 

2400 

12000 

2400 

6000 

1200 

8 000 

1600 

55 600 

17 600 

Survey 

Full 

Half-Ring 



7 200 

14400 

8 000 

1600 

20000 

4000 

24000 

4 800 

24000 

4 800 

12000 

2400 

16000 

3 200 

111200 

35 200 

Total 

72000 

72000 

118 800 

41400 

73 800 

84600 

84 600 

52200 

72000 

671400 


FFP8 CMB MC Maps 

Detectors 

Interval 

Cut 

Frequency 

Total 

30 

44 

70 

100 

143 

217 

353 

545 

857 

Channel 

Mission 

Full 

10000 

10000 

10000 

10000 

10000 

10000 

10000 

10000 

10000 

90000 

Half-Mission 

Full 

20000 

20000 

20000 

20000 

20000 

20000 

20000 

20000 

20000 

180000 

Detset 

Mission 

Full 



30000 

20000 

50000 

40000 

20000 

20000 

10000 

190000 

Total 

30000 

30000 

60000 

50000 

80000 

70000 

50000 

50000 

40000 

460000 


FFP8.1 CMB MC Maps 

Detectors 

Interval 

Cut 

Frequency 

Total 

30 

44 

70 

100 

143 

217 

353 

545 

857 

Channel 

Mission 

Full 

1000 

1000 

1000 

1000 

1000 

1000 

1000 

1000 

1000 

9 000 

Half-Mission 

Full 

2000 

2000 

2000 

2000 

2000 

2000 

2000 

2000 

2000 

18 000 

Total 

3 000 

3 000 

3 000 

3 000 

3 000 

3 000 

3 000 

3 000 

3 000 

27 000 
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Table 2. Values of tensor (r) and non-Gaussianity (/nl) param¬ 
eters for the FFP8 baseline and FFP8a-d blind challenge cases 


Data set 

r 

/nl 

FFP8. 

0 

0 

FFP8a . 

0.222 

-4.311 

FFP8b . 

0.046 

8.590 

FFP8c . 

0.088 

7.140 

FFP8d . 

0.153 

-2.181 


While the full set of maps can be generated for the fidu¬ 
cial cases, for the 10^-realization MC sets this would result in 
some 10^ maps and require about 6 PB of storage. Instead, there¬ 
fore, the number of realizations generated for each type of map 
is chosen to balance the improved statistics it supports against 
the computational cost of its generation and storage. In practice, 
we generate 10^ MC maps in total, including the full 10"^ re¬ 
alizations for all channel/mission/full maps for both noise and 
CMB. As detailed in Table 1, the remaining noise MCs sample 
broadly across all data cuts, while the additional CMB MCs are 
focused on the channel/half-mis sion/full maps and the subset of 
the detset/mission/full maps required by the Commander compo¬ 
nent separation code (Planck Collaboration X 2015). 

3. Inputs 

The inputs to an FFP simulation consist of the characterization 
of the mission (the satellite pointing and each detector’s data 
fiags), the intrument (the focal plane layout and the noise statis¬ 
tics, beam, and bandpass of each of the detectors), and the syn¬ 
thetic sky to be observed. 

3.1. Mission and instrument characteristics 

The goal of FFP8 is to simulate the Planck mission as accu¬ 
rately as possible; however, there are a number of known sys¬ 
tematic effects that are not included, either because they are re¬ 
moved in the pre-processing of the time-ordered data (TOD), or 
because they are insufficiently well-characterized to simulate re¬ 
liably, or because their inclusion (simulation and removal) would 
be too computationally expensive. These systematic effects are 
discussed in detail in Planck Collaboration II (2015) and Planck 
Collaboration VII (2015) and include: 

“ cosmic ray glitches (HFI); 

“ spurious spectral lines from the 4-K cooler electronics (HFI); 

- nonlinearity in the analogue-to-digital converter (HFI); 

- imperfect reconstruction of the focal plane geometry. 

Note that if the residuals from the treatment of any of these ef¬ 
fects could be mapped in isolation, then maps of such systemat- 
ics could simply be added to the existing FFP8 maps to improve 
their correspondence to the real data. 

3.1.1. Pointing and flags 

The FFP8 detector pointing is calculated by interpolating the 
satellite attitude to the detector sample times and by applying 
a fixed rotation from the satellite frame into the detector frame. 
The fixed rotations are determined by the measured focal plane 
geometry (Planck Collaboration IV 2015; Planck Collaboration 
VII 2015), while the satellite attitude is described in the Planck 


attitude history files (AHF). The FFP pointing expansion repro¬ 
duces the DPC pointing to sub-arcsecond accuracy, except for 
three short and isolated instances during Surveys 6-8 where the 
LFI sampling frequency was out of specification. Pixelization 
of the information causes the pointing error to be quantized to 
either zero (majority of cases) or the distance between pixel cen¬ 
tres (3.'4 and 1.'7 for LFI and HFI, respectively). Since we need a 
single reconstruction that will serve both instruments efficiently 
in a massively parallel environment, we use the pointing pro¬ 
vided by the Time Ordered Astrophysics Scalable Tools (TOAST) 
package. 

In Fig. 1 we compare LFI DPC and TOAST “hit maps,” which 
show the number of times a detector beam centre has fallen in a 
given pixel at the mid-point of a data sample, for the specified 
time period. There is a low level, diffuse, difference during the 
first three surveys, caused by slight differences in the way the 
pointing correction is interpolated in the two pipelines. The net 
effect is that some samples near pixel boundaries are assigned 
to different pixels in the different pipelines. Due to pixels being 
vastly oversampled, this affects about 25 % of the sky pixels in 
the first three surveys, while the fraction of discrepant pixels af¬ 
ter Survey 3 is about 0.5 %. In the affected pixels, the relative 
difference in hit count is 0.2 + 0.1 % throughout. For Survey 1, 
the 70 GHz median hit count is 812 and 68 % of the hit counts 
fall between 604 and 1 278. 

Except for Survey 7, the channel integration time differs by 
about a second, with the flight data maps always containing 
slightly more data. The larger difference during LFI Survey 7 
occurs where TOAST rejected 2.5 min of data due to their irregu¬ 
lar sampling rate. All of these differences are entirely negligible 
compared to the 14 Ms of overall integration time in each survey. 
For Survey 7 the relative difference in hits in differing pixels is 
slightly elevated to 0.3 %, with an asymmetric standard devia¬ 
tion of 0.6 %. 

The survey-by-survey HFI comparison in Fig. 2 does not 
have a diffuse pattern matching Fig. 1 that ends during Survey 
3, indicating a better agreement in the pointing correction inter¬ 
polation between TOAST and HFI DPC. The more pronounced 
ecliptic pole caustics are caused by a truncation error in the TOD 
time stamps that were used as a basis of pointing expansion. In 
addition, the HFFs smaller pixel size and higher sampling fre¬ 
quency increases the likelihood that small numerical differences 
can result in samples being assigned across a pixel boundary. On 
average, about 52 % of the pixels are affected and exhibit hit dif¬ 
ferences at the level of 0.5 + 0.3 %. For 100 GHz and Survey 1, 
the median hit count is 224 and 68 % of the hit counts fall be¬ 
tween 164 and 361. The channel integration time between the 
two pipelines in the HFI case shows negligible differences, of 
the order of a few tenths of a second for all surveys. 

3.1.2. Noise 

We require simulated noise realizations that are representative 
of the noise in the flight data, including variations in the noise 
power spectral density (PSD) of each detector over time. To ob¬ 
tain these we developed a noise estimation pipeline complemen¬ 
tary to those of the DPCs. The goal of DPC noise estimation is 
to monitor instrument health and to derive optimal noise weight¬ 
ing, whereas our estimation is optimized to feed into noise sim¬ 
ulation. Key features are the use of full mission maps for sig¬ 
nal subtraction, long (about 24 hour) realization length, and the 
use of autocorrelation functions in place of Fourier transforms 
to handle flagged and masked data (HFI). 
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Survey 1 Survey 2 







-4 Ahits 4 

Fig. 1. Survey hit map differences between FFP8 (using TOAST 
pointing) and LFI at 70 GHz. The diffuse error pattern during 
Surveys 1-3 was caused by differences in pointing interpolation. 
The matching HFI comparison is presented in Fig. 2. The FFP 
pipeline rejected 2.5 min of integration time on Survey 7, leaving 
a small deficit of hits on a single scanning ring that is not visible 
in this projection. A hit is assigned to a pixel where the detector 
beam centre falls midway through a sample integration period. 

The PSD for each detector and each roughly 45 min pointing 
period was estimated as follows. 

1. Read in a 15% mask for the Galaxy and point sources. 

2. Read the calibrated and dipole-subtracted time-ordered data 
and quality flags for the interval. 

3. Expand detector pointing for the interval. 

4. Scan and subtract a signal estimated from the full mission, 
full frequency map using the detector pointing weights. 

5. Subtract the mean of the interval. 

6. Then either 

- Fourier transform the signal-subtracted TOD for a noise 
PSD estimate, correcting for missing samples by divid¬ 
ing by the unflagged fraction (LFI), 

or 

- estimate the autocovariance function from the unflagged 
and unmasked samples and Fourier transform the auto¬ 
covariance function into a noise PSD (HFI). 

7. Bin the noise PSD estimate into 300 logarithmically spaced 
bins. 

8. Average the noise estimates into approximately daily noise 
spectra. 


Survey 1 Survey 2 



1 ■ 

-4 Ahits 4 


Fig. 2. Survey hit map differences between FFP8 and HFI at 
100 GHz. The diffuse pattern is caused by time-stamp truncation 
and is not the same as the pointing interpolation difference in the 
LFI comparison of Surveys 1-3 in Fig. 1. 


9. Fit an analytic 1 // model to the measured spectra. 

We show the results of running our pipeline on a signal-l-noise 
simulation in Fig. 3. 



Frequency [Hz] 

Fig. 3. Example comparison of the input and recovered noise 
PSDs for one pointing period of a signal-hnoise simulation of 
HFI bolometer 100-la. The shaded regions around the input 
model reflect the asymmetric realization scatter of the estimated 
PSDs at 68 % and 95 % confidence intervals in each of the 707 
logarithmically-placed frequency bins. 

As a part of this process we observed significant low fre¬ 
quency (~ 1 mHz) correlated noise in the flight data between 
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Fig. 4. BB spectra of 70 GHz {left) and 100 GHz {right) half-ring half difference noise maps, for both FFP8 simulations {green) 
and flight data {black). Differences at the low-^ end are caused by sample variance. These are pseudo-spectra, computed on 15% 
of the sky with the Galactic plane and point sources masked. Top: Linear horizontal axis to show small-scale behaviour. Bottom: 
Logarithmic horizontal axis to show large-scale behaviour. 


polarized HFI bolometers that share a horn,^ designated a and 
b. These pairs of detectors form a polarization-sensitive unit, as 
their relative difference is sensitive to linear polarization of the 
incident radiation, and the existence of a correlated noise com¬ 
ponent between them is a known feature of the HFI intrument 
Planck Collaboration VI (2014). When this correlated noise was 
not included in the simulations, it left a clear signature in the 
mismatch between our noise simulations and the half ring differ¬ 
ence maps. Specifically, ajb correlated noise adds coherently in 
the {a -r b) combination that is sensitive to temperature modula¬ 
tion of the signal, whereas it cancels in the {a - b) combination 
that detects polarization modulation of the signal. Without simu¬ 
lating the correlated noise component, our simulated noise maps 
had too little large-scale noise in the TT spectrum, but excess 
noise in the EE and BB spectra. We estimated the correlated 
noise spectrum by differencing the {a b) and {a - b) noise esti- 


^ Each of the polarized horns on the Planck focal plane has two or¬ 
thogonally polarized detectors. Due to the scanning strategy, it takes two 
polarized horns, i.e., a quadruplet of detectors, to solve for polarization 
away from the regions close to the ecliptic poles where the scanning 
rings intersect. 


mates, and added the correlated component into our noise simu¬ 
lation. FFP8 is the first simulation to include realistic ajb noise 
correlations with the appropriate red PSD. 

We also found that our physically motivated HFI noise mod¬ 
els contained a mismatch with the measured noise PSDs at the 
very highest frequencies (> 50 Hz). Replacing the models with 
actual binned noise PSDs in the highest-frequency bins gave bet¬ 
ter performance in the half-ring map test. Examples of the corre¬ 
lated noise spectrum and the noise model fits to the autospectra 
can be found in Planck Collaboration VII (2015). 

To validate the results we considered the angular power spec¬ 
tra of both simulated and flight data half-ring map differences 
(i.e., splitting each pointing period into two halves, making maps 
from all of the first halves and second halves separately, and tak¬ 
ing the difference of the maps), chosen because they are virtually 
free of all signal residuals while still supporting most of the same 
noise power as the full-ring data. Fig. 4 illustrates the excellent 
agreement between the simulated and real data half-ring map BB 
spectra for 70 and 100 GHz. 
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3.1.3. Beams 

The simulations use the so-called scanning beams (e.g., Planck 
Collaboration VI 2014), which give the point-spread function 
of for a given detector including all temporal data process¬ 
ing effects: sample integration, demodulation, ADC nonlinearity 
residuals, bolometric time constant residuals, etc. In the absence 
of significant residuals (LFI), the scanning beams may be esti¬ 
mated from the optical beams by “smearing” them in the scan¬ 
ning direction to match the finite integration time for each instru¬ 
ment sample. Where there are unknown residuals in the time¬ 
lines (HFI), the scanning beam must be measured directly from 
observations of strong point-like sources, namely planets. If the 
residuals are present but understood, it is possible to simulate the 
beam measurement and predict the scanning beam shape starting 
from the optical beam. 

For FFP8, the scanning beams are expanded in terms of their 
spherical harmonic coefficients, with the order of the expan¬ 

sion (maximum I and m considered) representing a trade-off be¬ 
tween the accuracy of the representation and the computational 
cost of its convolution. The LFI horns have larger beams with 
larger sidelobes (due to their location on the outside of the focal 
plane), and we treat them as full An beams divided into main (up 
to 1?9, 1?3, and 0?9 for 30, 44, and 70 GHz, respectively), inter¬ 
mediate (up to 5°), and sidelobe (above 5°) components (Planck 
Collaboration IV 2015). This division allows us to tune the ex¬ 
pansion orders of the three components separately. HFI horns 
are limited to the main beam component, measured out to 100' 
(Planck Collaboration VII2015). Since detector beams are char¬ 
acterized independently, the simulations naturally include differ¬ 
ential beam and pointing systematics. 


3.1.4. Bandpasses 

Both the LFI and HFI detector bandpasses are based on ground 
measurements (see Zonca et al. 2009 and Planck Collaboration 
IX 2014, respectively), although flight data processing for LFI 
now uses in-flight top-hat approximations rather than the ground 
measurements that were found to contain systematic errors. 
Differences in the bandpasses of detectors nominally at the same 
frequency (the so-called bandpass mismatch) generate spurious 
signals in the maps, since each detector is seeing a slightly differ¬ 
ent sky while the mapmaking algorithms assume that the signal 
in a pixel is the same for all detectors. To quantify the effect 
of these residuals, in FFP8 we generate detector timelines from 
foreground maps in two ways, one that incorporates the individ¬ 
ual detector bandpasses, the other using an average bandpass for 
all the detectors at a given frequency. 

This effect of the bandpass mismatch can be roughly mea¬ 
sured from either flight or simulated data using so-called spuri¬ 
ous component mapmaking, which provides noisy all-sky esti¬ 
mates of the observed sky differences (the spurious maps), ex¬ 
cluding polarization, between individual detectors and the fre¬ 
quency average. We compare the amount of simulated bandpass 
mismatch to flight data in Figs. 5 (LFI) and 6 (HFI). The spurious 
component approach is detailed in Appendix A. Mismatch be¬ 
tween FFP8 and flight data is driven by inaccurate bandpass de¬ 
scription (LFI) and incomplete line emission simulation (HFI). 
The noisy pixels that align with the Planck scanning rings in 
the HFI maps are regions where the spurious map solution is 
degenerate with polarization due to insufficient observation ori¬ 
entations. 


3.2. The foreground sky 

The quality of the reconstructed CMB signal depends on the ac¬ 
curacy of the separation of the astrophysical foreground com¬ 
ponents from the Planck data, which requires building simu¬ 
lations as close to the complexity of the real sky as possible 
over a large frequency range. To this end, a complete model 
of multi-component sky emission, called the Planck sky model 
(PSM), has been developed. FFP8 simulations make use of ver¬ 
sion 1.9.0 of the PSM, a major update of the published ver¬ 
sion 1.7.8 (Delabrouille et al. 2013). The foreground compo¬ 
nents generated by the PSM are diffuse Galactic emission from 
thermal dust, spinning dust, relativistic electrons (synchrotron 
emission), free-free radiation, and CO (the /=1^0, 7=2^1, and 
J=3^2 lines at 115.27, 230.54, and 345.80 GHz, respectively), 
plus the cosmic infrared background (CIB), emission from radio 
sources, and the thermal and kinetic Sunyaev-Zeldovich (SZ) ef¬ 
fects. 

3.2.1. Thermal dust 

The thermal dust emission is modelled using single-frequency 
template maps of the intensity and polarization, together with a 
pixel-dependent emission law. 

For FFP8 the thermal dust emission templates are derived 
from the Planck 353 GHz observations. This update of the origi¬ 
nal PSM dust model, described in Delabrouille et al. (2013) and 
based on Finkbeiner et al. (1999), is necessary to provide a bet¬ 
ter match to the emission observed by Planck. While one option 
would be simply to use the dust opacity map obtained in Planck 
Collaboration XI (2014), this map still suffers from significant 
contamination by CIB anisotropies and infrared point sources. 
Using it as a 353 GHz dust template in simulations would re¬ 
sult in an excess of small scale power (from CIB and infrared 
sources) scaling exactly as thermal dust across frequencies. The 
resulting component represents correctly neither dust alone (be¬ 
cause of an excess of small scale power) nor the sum of dust and 
infrared sources (because the frequency scaling of the CIB and 
infrared sources is wrong). For simulation purposes, the main 
objective is not to have an exact map of the dust, but instead a 
map that has the right statistical properties. Hence we produce a 
template dust map at 353 GHz by removing that fraction of the 
small-scale power that is due to CIB emission, infra-red sources, 
CMB, and noise. 

The emission of strong point sources is first subtracted from 
all Planck 2015 maps and from the IRAS 100 yum map following 
Miville-Deschenes & Lagache (2005). We mask point sources 
detected in Planck and IRAS maps with an algorithm designed 
specifically to avoid false detections in regions of dust emis¬ 
sion where localized compact regions of diffuse dust emission 
can easily be mistaken for point sources. This algorithm detects 
sources against a background level computed directly on the 
maps in the immediate vicinity of each pixel. Specifically, for 
each pixel p, we compute the standard deviation (T{p) of the ob¬ 
served map in a neighbourhood N{p) consisting of a small ring 
centred on the source, between the two radii ri{p) and r 2 (p). We 
consider as part of point sources those objects for which the flux 
is above a threshold of 3 to 5 times (T{p). The parameters of 
the filter for each frequency map - the full-width half-maximum 
(FWHM) radii defining the disc, and the source detection thresh¬ 
old - are optimized by trial and error. We then generate a mask 
from the union of small discs of radius 1 FWHM centred on all 
pixels above the detection threshold at any frequency. The gaps 
created by masking the point sources are filled by extrapolation 
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Fig. 5. Pairs of measured bandpass mismatch LFI maps, comparing FFP8 simulations (upper) and 2015 Planck (lower) data. Our 
method of estimating the total bandpass mismatch between each detector and a frequency average is detailed in Appendix A. 


of the neighbouring regions using a minimum-curvature spline 
gap filler. This procedure avoids most false detections of point 
sources in filaments of dust emission, unless there is a compact 
source significantly brighter than the local emission of the fila¬ 
ment. The drawback of this approach is that extended compact 
sources such as nearby galaxies are not masked. 

The dust intensity template is then corrected for CIB, CMB, 
and noise fluctuations using the multi-dimensional Generalized 
Needlet Internal Linear Combination (GNILC) of Remazeilles 
et al. (201 1). Starting from the 10 de-sourced maps (nine Planck 
maps and the IRAS lOOyum map), we model the total emission 
at a given frequency v as: 

( 1 ) 


where yy is the observed map, ^-y the map of foreground emission 
(mostly Galactic emission), and Xy the sum of CMB, CIB, and 
noise. The total 10 x 10 covariance matrix Ry of the maps, in 
any localized region in pixel and/or harmonic space, is the sum 
of the foreground term Rs and a contamination term R^ due to 
CMB, CIB and noise, 

Ry = Rs + Rx- (2) 

We perform a needlet decomposition of each map as in 
Remazeilles et al. (2011), except that the spectral windows cor¬ 
respond to Gaussian bandpass filters in harmonic space. We use 
the best fit CMB Q from Planck Collaboration XV (2014), the 
CIB power spectrum at each frequency and cross-correlation co¬ 
efficients from Planck Collaboration XXX (2014), and the noise 
level from Planck half-ring difference maps to obtain a model Rx 
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Fig. 6. Pairs of measured bandpass mismatch HFI maps, comparing FFP8 simulations (upper) and 2015 Planck (lower) data. Our 
method of estimating the total bandpass mismatch between each detector and a frequency average is detailed in Appendix A. 
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of the covariance Rx of CMB, CIB, and noise in a set of needlet 
domains localized on the sky. We then follow the procedure de¬ 
scribed in Remazeilles et al. (2011) to reconstruct maps of at 
up to 5' resolution. For each needlet domain, the covariance of 
the needlet power is whitened locally as 

Ry = (3) 


The matrix Ry is diagonalized and the Akaike information crite¬ 
rion (AIC; Akaike 1974) is used to filter out the CMB, CIB, and 
noise contamination Xy in dimensions of the observation space 
where the total covariance is compatible with ^-y = 0. The final 
map of foreground emission ^-y at each needlet scale is recon¬ 
structed by coadding the needlets for which there is significant 
additional emission from Galactic foregrounds according to the 
AIC criterion. However, the original source-subtracted 353 GHz 
Planck intensity map is used instead of the processed map in 
the Galactic plane, the LMC, the SMC, regions of very strong 
localized emission, and the extended galaxies M31, M33, M81, 
M82, M87, MlOl, NGC 55, NGC 300, and NGC 2403. The rea¬ 
son for this is that the multi-dimensional GNILC process is not 
perfectly local in pixel space, so that bright local structures such 
as the Galactic plane or nearby galaxies are analysed together 
with regions of significantly lower emission. This results in ex¬ 
cessive filtering of some small-scale structure in the bright ob¬ 
jects themselves, and insufficient filtering in the pixels located 
in their immediate vicinity. While the latter is not much of a 
problem (giving a non-crucial loss of efficiency of the multi¬ 
dimensional GNILC filtering), the former is not acceptable since 
it impacts the brightest objects in the sky. As CIB, CMB, and 
noise contributions are completely subdominant in such regions 
of very bright emission, it is better to use the original map than 
the processed one for the generation of the dust template in these 
regions. 

Next, a mask is applied to remove from the 353 GHz dust 
map residuals of the positive SZ effect towards bright galaxy 
clusters. The mask consists of discs centred on CLASS 1 (i.e., 
the most reliable) SZ candidates in the PSZl catalogue (Planck 
Collaboration XXIX 2014), with radius 5 x OsqqJ Gaps gen¬ 
erated in the map by this masking are filled in again by 
a minimum-curvature spline “inpainting”. Finally, a colour- 
correction factor is applied to obtain a map of brightness emis¬ 
sion at 353 GHz rather than the map integrated in the Planck 
353 GHz channel. The 353 GHz dust map obtained in this 
way, which is subsequently used as a dust template, is illus¬ 
trated in Fig. 7, where it is also compared with the model of 
Finkbeiner et al. (1999) and with the map of optical depth pub¬ 
lished in Planck Collaboration XI (2014). 

The pixel-dependent emission law is derived from Planck 
Collaboration XI (2014), fitting Planck 353, 545, and 857 GHz 
and IRAS 100 jam data with a modified blackbody. The emission 
of the strongest point sources is subtracted from all maps before 
performing the fits, with gaps filled by the same method used for 
the dust intensity template. 

For each pixel p, the spectral energy distribution is approxi¬ 
mated with a modified blackbody, i.e.. 


Ivip) = I 


353 GHz 


(/^)( 


353 GHz 



BvjTdnstip)) \ 
B353Gliz(Tdusl(p)) I 


(4) 


where the maps of the dust temperature T^ust and spectral index 
ySpiR are those presented in Planck Collaboration XI (2014). The 


^ 6500 is the observed radius Rsoo/Da(z), where Da(z) is the angular 
diameter distance of the cluster. 



-0.1 -0.03 0.03 0.1 

MJy/sr 


Fig. 7 . 12?5 X 12?5 gnomonic projection of dust intensity tem¬ 
plates centred at high latitude (l,b) = (90°,-80°) at a map 
resolution of 5'. Top: new PSM dust from multi-dimensional 
GNILC (Remazeilles et al. 2011). Middle: original PSM dust 
(Delabrouille et al. 2013) based on the model of Finkbeiner et al. 
(1999). Bottom: Planck 2013 dust model (Planck Collaboration 
XI 2014). Comparing the top and bottom panels, the improve¬ 
ment in reducing the contamination of the dust by small-scale 
CIB anisotropies at high Galactic latitude is clearly visible. 


10 









Planck Collaboration: Full Focal Plane Simulations 


sky map is produced by applying the multi-dimensional 

GNILC of Remazeilles et al. (201 1). 

Maps of the Qjl and Ujl ratios between Stokes parameters 
were computed from the best versions of the Planck maps avail¬ 
able at the time of generation of the simulations in June 2014, 
smoothed to 30' resolution to reduce noise. At 353 GHz, the ob¬ 
served maps are used for polarization, while the intensity map is 
obtained from multi-dimensional GNILC processing. All maps 
are corrected for the dust bandpass mismatch leakage. Structure 
on small scales is added in order to have no break in the power 
spectra of dust polarization. These maps determine the polariza¬ 
tion orientation and fraction. For each frequency, the Stokes Q 
and U maps of the model are then obtained by multiplying the 
Qjl and UII maps by the total intensity map at full resolution. 
The polarization orientation does not depend on frequency, but 
the difference between and ySpiR introduces a spectral depen¬ 
dence of the polarization fraction, which varies over the sky. 

To extrapolate the model to frequencies below 353 GHz we 
follow Planck Collaboration Int. XVII (2014), using a spectral 
index map y^mm slightly smaller than ySpiR for the total dust in¬ 
tensity, but not for the polarized emission. The difference be¬ 
tween ySmm and y^FiR IS a quadratic function of r^ust. which has 
been fitted comparing the Planck maps at 143 and 353 GHz, after 
subtraction of the CMB anisotropies. It introduces a small spec¬ 
tral difference between dust intensity and polarization that varies 
over the sky, similar to that reported in Planck Collaboration Int. 
XXII (2015). 


3.2.2. Spinning dust 

Several observations indicate an excess of microwave emission 
in the 10—100 GHz range (Kogut et al. 1996; Leitch et al. 1997; 
de Oliveira-Costa et al. 1999; Watson et al. 2005). This emission 
has been proposed to be the dipole radiation from rapidly spin¬ 
ning nano-scale dust grains (Draine & Lazarian 1998). Detailed 
analyses of the data from various experiments show detections 
of a component that is consistent with the spinning dust model 
(Davies et al. 2006; Miville-Deschenes et al. 2008; Kogut et al. 
2011; Planck Collaboration XX 2011; Planck Collaboration Int. 
XV 2014; Planck Collaboration X 2015). 

The spinning dust map used for FFP 8 simulations is a simple 
realization of the spinning dust model, post-processed to remove 
negative values occurring in a few pixels because of the genera¬ 
tion of small-scale fluctuations on top of the spinning dust tem¬ 
plate extracted from WMAP data by Miville-Deschenes et al. 
(2007). The model uses a single template map, illustrated in Fig. 
8 , which is scaled by a single emission law parameterized fol¬ 
lowing the model of Draine & Lazarian (1998). 


3.2.3. Synchrotron 

At low frequencies (below about 80 GHz) the dominant contam¬ 
inant of the polarized CMB signal is the Galactic synchrotron 
emission that arises from the acceleration of relativistic cosmic 
rays in the Galactic magnetic field. FFP 8 synchrotron emission 
is modelled on the basis of the template emission map observed 
at 408 MHz by Haslam et al. (1982). This template synchrotron 
map is extrapolated in frequency using a spectral index map cor¬ 
responding to a simple power law. 



0 MJy/sr 60 

Fig. 8. Template of spinning dust intensity. 


3.2.4. Free-free 

Electron-ion interactions in the ionized phase of the ISM pro¬ 
duce emission that is in general fainter than both the synchrotron 
and the thermal dust emission outside of the active star-forming 
regions in the Galactic plane. The free-free model uses a single 
template, which is scaled in frequency by a specific emission 
law. The free-free spectral index is a slowly varying function of 
frequency and depends only slightly on the local value of the 
electron temperature Te (Bennett et al. 1992, 2003; Dickinson 
et al. 2003). In FFP 8 , the free-free spectral dependence is mod¬ 
elled by assuming a constant electron temperature Te = 7000 K, 
as described in Dickinson et al. (2003). 

3.2.5. Cosmic infrared background 

The GIB model relies on the distribution of individual galax¬ 
ies in template maps based on the distribution of dark matter 
at a range of relevant redshifts. We assume the GIB galaxies 
can be grouped into three different populations (proto-spheroid, 
spiral, starburst) as described in Cai et al. (2013). Within each 
population, galaxies have the same SED, while the flux den¬ 
sity is randomly distributed according to redshift-dependent 
number counts obtained from JCMT/SCUBA-2 observations 
(Chen et al. 2013) and the Planck ERCSC (Negrello et al. 
2013), as well as observations cited in Cai et al. (2013), most 
notably Herschelj^^ixt (Hermes; Bethermin et al. 2013) and 
AzTEC/ASTE (Combined ACES; Scott et al. 2012). 

We use the Class (Bias et al. 2011) software to generate 
dark matter maps at 17 different redshifts between 1 and 5.5. 
Since the galaxy distribution does not exactly follow the dark 
matter distribution, we modify the a^^ coefficients of dark mat¬ 
ter anisotropies given by Class at redshift Zi such that 

atmjiZi) = X bjizi) x |l + ||| j, (5) 

where bj(z) is the population-dependent component of the bias 
of the population 7 , and is identical for spirals and starbursts. 
(1 -r (I) ) is the scale-dependent part of the bias. Values of 
a = 1.44 and JS = 2 654 were obtained by fitting Planck spectra. 
Template maps generated from the a^mjizi) coefficients are then 
exponentiated to avoid negative pixels. 

Galaxies are randomly distributed with a probability of pres¬ 
ence proportional to the pixel values of the template maps. One 
map is generated for each population, at each redshift, and asso¬ 
ciated with a redshifted SED depending on the population. The 
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Fig. 9. Simulated and observed cosmic infrared background 
power spectra. Solid lines are from the PSM simulations; dashed 
lines are the Planck “extended halo model” which includes the 
correlated and Poisson contributions to the CIB power spec¬ 
trum from dusty galaxies. Data points are Planck observations, 
corrected for the contribution of radio sources, all from Planck 
Collaboration XXX (2014). 


emission of these maps (initially at a reference frequency) can 
be extrapolated to any frequency using the associated redshifted 
SED. By summing the emission of all maps, we can generate 
CIB maps at any frequency in the range of validity of our model. 

A comparison of CIB power spectra of the simulated maps 
and the maps obtained by Planck observations at HFI frequen¬ 
cies is shown in Fig. 9; agreement is excellent at all but the 
lowest two HFI frequencies where the CIB is subdominant. See 
Planck Collaboration XXX (2014) for further discussion. 

3.2.6. Radio sources 

Radio source modelling in the PSM was described in detail 
in Delabrouille et al. (2002). Here, we describe only the main 
changes since that pre-launch version. 

For “strong” radio sources (5 30 > 0.5 Jy), we use radio 
sources compiled by Gonzalez-Nuevo et al. (2008) at 0.84, 1.4, 
or 4.85 GHz. For sources observed at two of these frequencies, 
we extrapolate or interpolate to the third frequency assuming 
the spectral index estimated from two observed. For sources ob¬ 
served at only one frequency, we use differential source counts 
from de Zotti et al. (2010) to obtain the ratio of steep- to flat- 
spectrum sources in each interval of flux density considered. 
From this ratio, we assign spectral indices (randomly) to each 
source within each flux density interval. Fiducial Gaussian spec¬ 
tral index distributions as a function of spectral class are ob¬ 
tained from the literature (Massardi et al. 2011; Sadler et al. 
2008). These are then adjusted slightly until there is reason¬ 
able agreement between the PSM differential counts and the 
model counts predicted by de Zotti et al. (2010). The final 
means and standard deviations of these spectral index distribu¬ 
tions are - 0 .7+0.2 for steep-spectrum sources, -0.3± 0.2 for flat- 
spectrum radio quasars (FSRQ), and -0.1 + 0.2 for BF Lacertae 
objects (BF Lacs). 

For “faint” radio sources (5 30 < 0.5 Jy), the pre-launch PSM 
showed a deficit of sources resulting from inhomogeneities in 
surveys at different depths. We address this issue by construct¬ 


ing a simulated catalogue of sources from the model differen¬ 
tial counts of de Zotti et al. (2010) at 1.4 GHz. We replace the 
simulated sources by the observed ones, wherever possible. If, 
however, in any particular pixel, we have a shortfall of observed 
sources, we make up the deficit with the simulated sources. 
Every source in this new catalogue is given a model-derived 
spectral class. We thus assign a spectral index to each source 
based on the spectral class, and model the spectrum of each 
source using four power laws covering the ranges v < 5 GHz, 
5 < y < 20 GHz, 20 < v < 100 GHz, and v > 100 GHz. We also 
assume some steepening of the spectral index with frequency, 
with fiducial values of the steepening obtained from the litera¬ 
ture (Massardi et al. 2011; Sadler et al. 2008). The steepening 
of the spectral index assumed for 5-20 GHz is -0.3 for steep- 
spectrum sources and -0.3 for flat-spectrum sources (FRSQs 
and BL Lacs). For 20-100 GHz, we assume no further steep¬ 
ening for steep-spectrum sources, and steepening of -0.3 for 
FRSQs and BL Lacs. 

We combine the faint and strong radio source catalogues we 
constructed and compute the differential source counts on these 
sources between 0.005 Jy and 1 Jy (see Fig. 10). These show 
good agreement with the model differential counts predicted by 
de Zotti et al. (2010), particularly at the faintest end. Finally we 
also model the polarization of these radio sources using the mea¬ 
sured polarization fractions from Ricci et al. (2004), compris¬ 
ing 71 steep- and 126 flat-spectrum sources; for each simulated 
source we draw a polarization fraction at random from the list of 
real sources of the same spectral type. 

3.2.7. SZ clusters 

SZ clusters are simulated following the model of Delabrouille, 
Melin, and Bartlett (DMB) as implemented in the PSM 
(Delabrouille et al. 2002). A catalogue of halos is drawn from 
a Poisson distribution of the mass function from Tinker et al. 
(2008) with a limiting mass of M 5 oo,tme > 2 x lO^^M©.^ We 
use the pressure profile from Arnaud et al. (2010) to model the 
thermal SZ emission of each halo given its redshift and mass, 
M 5 oo,x- In order to match the observed counts in Planck, we set 
4^500,x/ 4 / 500 ,true = 0-65. Fixing this ratio to unity would lead to 
twice as many clusters simulated as appear in the actual data, 
as shown in Planck Collaboration XX (2014). We determine the 
cluster temperature from the relation of Arnaud et al. 

(2005), and assume that the profiles are isothermal. These steps 
allow us to compute the first-order thermal relativistic correc¬ 
tion (Itoh et al. 1998) and the kinetic SZ effect for each cluster, 
both of which are included in the simulation. Finally, we inject 
MCXC (Piffaretti et al. 201 1) and MaxBCG (Koester et al. 2007) 
clusters following the same model, and remove from the simu¬ 
lation corresponding clusters in each redshift and mass range. 
Hence the SZ simulation features the majority of known X-ray 
and optical clusters, and is fully consistent with X-ray scaling 
laws and observed Planck SZ counts. 

3.2.8. Total foreground sky 

The sky model is simulated at a resolution common to all com¬ 
ponents by smoothing the maps with an ideal Gaussian beam 


^ M 500 is the total cluster mass within a sphere of radius R 5 Q 0 , defined 
as the radius within which the mean mass over-density of the cluster 
is 500 times the cosmic critical density at the redshift of the cluster. 
Msm = (4;r/3)R5oo3[500pc(z)], with pdz) = 3H^{z)K^7tG), where H{z) 
is the Hubble parameter with present-day value Hq. 
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log(S[Jy]) 



log(S[Jy]) 


Fig, 10. Counts of the faint and the strong radio sources at 
30 GHz (top) and 100 GHz (bottom). The differential counts of 
the faint and the strong radio sources between 0.005 Jy and 1 Jy 
are shown by vertical bars. 


of FWHM of 4'. The HEALPix^ pixelization in Galactic co¬ 
ordinates is used for all components, with A^side = 2048 and 
^max = 6000. Sky emission maps are generated by numerically 
band-integrating the sky model maps (emission law of each com¬ 
ponent, in each pixel) over the frequency bands both of each 
detector in the focal plane and - using an average over the detec¬ 
tors at a given frequency - of each channel. The band-integrated 
maps are essentially observations of the model sky simulated 
by an ideal noiseless instrument with ideal Gaussian beams of 
FWHM equal to the resolution of the model sky. 


3.3. The CMB sky 

The CMB sky is simulated in three distinct components, namely 
lensed scalar, tensor, and non-Gaussian complement. The total 
CMB sky is then the weighted sum with weights 1, V^, and /nl, 
respectively. For FFP8, all CMB sky components are produced 
as spherical harmonic representations of the /, Q, and U skies. 

The FFP8 CMB sky is derived from our best estimate of 
the cosmological parameters available at the time of its genera¬ 
tion, namely those from the first Planck data release (PR 1-2013; 
Planck Collaboration 12014), augmented with a judicious choice 
of reionization parameter r, as listed in Table 3. 

3.3.1. The scalar CMB sky 

The scalar component of the CMB sky is generated includ¬ 
ing lensing, Rayleigh scattering, and Doppler boosting effects. 
Using the camb code (Lewis et al. 2000), we first calculate fidu¬ 
cial unlensed CMB power spectra the lensing 

potential power spectrum and the cross-correlations 

and . We then generate Gaussian T, E, and 0 multipoles 
with the appropriate covariances and cross-correlations using a 
Cholesky decomposition and three streams of random Gaussian 
phases. These fields are simulated up to ^max = 5120. In addi¬ 
tion, we add a dipole component to 0 to account for the Doppler 
aberration due to our motion with respect to the CMB (Challinor 
& van Leeuwen 2002): 


(pp(h) = fin. 


( 6 ) 


where fi is in the direction of Galactic coordinates (/,/?) = 
(264?4,48?4) with amplitude \fi\ = 0.00123, corresponding to 
a velocity of 369 km s“^ (Hinshaw et al. 2009). 

We we compute the effect of gravitational lensing on the tem¬ 
perature and polarization fields, using an algorithm similar to 
LensPix (Lewis 2005). We use a fast spherical harmonic trans¬ 
form to compute the temperature, polarization, and deflection 
fields: 

[de + id^p\ = ^ - V^(^ + 1) 0^m- (J) 

The unlensed CMB fields T, Q, and U are evaluated on an 
equicylindrical pixelization (ECP) grid with Nq = 32768 and 
N^p = 65 536, while the deflection field is evaluated on a 
HEALPix A^side = 2048 grid. We then calculate the “lensed po¬ 
sitions” h' for each A^side = 2048 HEALPix pixel. In the fiat-sky 
limit these are given schematically as w' = w -r d(ji)\ in practice 
on the curved sky we parallel transport h a distance \d{h)\ along 
the great circle in the direction d{h). We then interpolate T, Q, 
U at the lensed positions using 2-D cubic Lagrange interpolation 
on the ECP grid. 

After lensing, we incorporate the frequency-dependent 
Doppler modulation effect (Planck Collaboration XXVII 
2014), multiplying the temperature and polarization maps by 
[l+byfi'h], with by factors given in Table 4. These boost 
factors incorporate bandpass corrections calculated given the 
frequency-dependence of the Doppler effect. We then evaluate 
lensed, Doppler boosted Eim, and up to ^max = 4096 
with a harmonic transform of the Aside = 2 048 HEALPix map of 
these interpolated T, Q, and U values. 

We next add frequency-dependent Rayleigh scattering ef¬ 
fects, using the “R” matrices defined in Eq. (7) of Lewis (2013). 

^ http://healpix.sourceforge.net 
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Table 3. Cosmological parameter values from the PRl-2013 and PR2-2015 data releases, used in FFP 8 and FFP8.1, respectively 


Parameter 

Symbol 

PRl-2013 (FFP8) 

PR2-2015 (FFP8.1) 

Baryon density. 

oj^ = 

0.0222 

0.0223 

Cold dark matter density. 

Cx)c - klch^ 

0.1203 

0.1184 

Neutrino energy density. 

CjL>y = 

0.00064 

0.00065 

Dark energy density. 

Da 

0.6823 

0.6931 

Hubble parameter, Hq = 100/zkms“^ Mpc"^ .... 

h 

0.6712 

0.6787 

Primordial curvature perturbation spectrum 
amplitude (at k - 0.05 Mpc"^). 

As 

2.09 X 10-9 

2.14 X 10-9 

spectral index . 

ns 

0.96 

0.97 

Thomson optical depth through reionization .... 

T 

0.065 

0.067 


Table 4. Boost factors 


Frequency 

[GHz] 

Boost factor 
br 

30. 

1.05 

44. 

1.1 

70. 

1.25 

100 . 

1.51 

143. 

1.96 

217. 

3.07 

353 . 

5.38 

545 . 

8.82 

857 . 

14.20 


We use effective central frequencies for the terms of 30, 44, 
70, 105, 148, 229, 372, 577, and 891 GHz for the 30-857 GHz 
channels, and effective central frequencies for the terms of 30, 
44, 70, 108, 151, 233, 378, 585, and 905 GHz. 

Finally we add a second-order temperature quadrupole 
following the equations in Kamionkowski & Knox (2003). 
Since the main Planck data processing removes the frequency- 
independent part (Planck Collaboration VIII 2015), we simulate 
only the residual frequency-dependent temperature quadrupole. 
After subtracting the frequency-independent part, the simulated 
quadrupole has frequency dependence oc (hy - l)/ 2 , which we 
calculate using the bandpass-integrated hy boost factors given in 
Table 4. Defining the z-axis to be in the direction of the solar 
dipole, at each frequency we add 


a 


T _ 
20 - 


4 (^ - 1 ) 
6 



( 8 ) 


3.3.2. The tensor CMB sky 

In addition to the scalar CMB simulations, we also generate a 
set of CMB skies containing primordial tensor modes. Using 
the fiducial cosmological parameters of Table 3, we calculate 

the tensor power spectra (yJTtensor^ tensor^ tensor 

ing camb with a primordial tensor-to-scalar power ratio r = 0.2 
at the pivot scale k = 0.05 Mpc“^ We then simulate Gaussian 
T, E, and 5-modes with these power spectra, and convert these 
to spherical harmonic representations of the corresponding /, Q 
and U maps. Note that the default r = 0.2 means that building 
t he FFP 8 a-d maps requires rescaling each CMB tensor map by 
Vr/0.2 for each of the values of r in Table 2. 


3.3.3. The non-Gaussian CMB sky 

We use a new algorithm to generate simulations of CMB 
temperature and polarization maps containing primordial non- 
Gaussianity. Non-Gaussian fields in general have a non¬ 
vanishing bispectrum contribution sourced by mode correla¬ 
tions. The bispectrum, the Fourier transform of the 3-point cor¬ 
relation function, can then be characterized as a function of three 
wavevectors, F(ki,k 2 ,k 3 ). Depending on the physical mecha¬ 
nism responsible for generating the non-Gaussian signal, it is 
possible to introduce broad classes of model that are categorized 
by the dependence of F on the type of triangle formed by the 
three momenta kf (see e.g., Babich et al. 2004). Here, we focus 
on non-Gaussianity of local type, where the bulk of the signal 
comes from squeezed triangle configurations, ki k 2 ^ k^. This 
is typically predicted by multi-field inflationary models (e.g., 
Moroi & Takahashi 2001; Enqvist & Sloth 2002; Lyth et al. 
2003). In this scenario, the primordial gravitational potential O 
is defined in terms of a Gaussian auxiliary held Ol, 

'I> = 'I>l+/nl('I>l-(®l>T (9) 

where higher-order terms are neglected (Verde et al. 2000; Wang 
& Kamionkowski 2000). Here, the non-Gaussian contribution to 
the primordial potential is parameterized by means of the scalar 
prefactor /nl- The spherical harmonic coefficients of the CMB 
anisotropies are then given by a linear equation, the line-of-sight 
integral (Liguori et al. 2003) 

^fm = ^ J' drr^^(„(r)af{r) 

= , (10) 

where we introduce the real space transfer function a^ir) for 
temperature (X = T) and polarization (X = E). We also make 
use of a compact matrix notation for the integral that, for all 
practical purposes, is evaluated numerically on a finite number 
of grid points. 

Our aim is to augment the Gaussian CMB simulations de¬ 
scribed above with non-Gaussian templates, defined by the sec¬ 
ond term in Eq. (9). Given the noiseless Gaussian map, we first 
compute a Wiener filter reconstruction of the primordial poten¬ 
tial Ol on the full sky. Let P be the covariance matrix of Ol on 
the radial grid (e.g., Liguori et al. 2003), 

/ A^ (k) 

dk (kri) 4 ikr2), (11) 
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where A^(^) is the primordial power spectrum and the are 
spherical Bessel functions. Further denoting C as the CMB co- 
variance matrix given by the fiducial power spectrum of the 
Gaussian simulated map i-, then 

(DWF(^) = (12) 

Since the Wiener filter is a biased estimate, we now supple¬ 
ment it with a fluctuation term to yield a Gaussian constrained 
realization (Bertschinger 1987). The covariance of the fiuctua- 
tion term is completely determined by the kernel of M, 

P^uc = p-PMtC“iMP. (13) 

We finally obtain an unbiased reconstruction of the linear 
gravitational potential for each simulation 5', 

OL(5)=<i)rw+(p‘’“")'^'^, ( 14 ) 

where g is a vector of univariate Gaussian random numbers. 
With the reconstruction in hand, it is now straightforward to ob¬ 
tain the non-Gaussian template 

^ng = M((Dl(^)2-((Dl)T (15) 

We evaluate all line-of-sight integrals (Eq. 10) on 70 radial 
shells, where quadrature weights and node positions are opti¬ 
mized as described in Eisner & Wandelt (2009). The resulting 
non-Gaussian map is exact in the sense that all higher-order cor¬ 
relation functions are correctly represented. 

3.3.4. The FFP8.1 CMB skies 

The FFP8 simulations are an integral part of the analyses used to 
derive PR2-2015, and so were necessarily generated prior to de¬ 
termining that release’s cosmological parameters. As such there 
is inevitably a mismatch between the FFP8 and the PR2-2015 
cosmologies, which we address in two ways. The quick-and- 
dirty fix is to determine a single rescaling factor that minimizes 
the difference between the PR 1-2013 and PR2-2015 TT power 
spectra and apply it to all of the FFP8 CMB maps; this number is 
determined to be 1.0134, and the rescaled maps have been used 
in several repeat analyses to confirm the robustness of various 
PR2-2015 results. 

More rigorously though, we also generate a second set of 
CMB realizations based on the PR2-2015 cosmology, dubbed 
FFP8.1, and perform our reanalyses using these in place of 
the FFP8 CMB skies in both the fiducial and MC realizations. 
Table 3 lists the cosmological parameters used for FFP8.1 while 
Table 1 enumerates the current status of the FFP8.1 CMB MCs. 

4. FFP8 pipelines 

The two Planck DPCs have distinct internal data formats and ac¬ 
cess protocols, including elements that are either impossible or 
impractical to port to supercomputing facilities. Moreover, data 
access can be a significant computational bottleneck, especially 
in an HPC environment, so aggressive 10 optimization is neces¬ 
sary for the most compute-intensive operations. To address these 
issues, we use the Planck exchange format^^ (EE), developed for 

The file format used in DPC exchanges is similar to the official 
released timelines (Planck Collaboration ES 2015) except that detec¬ 
tor pointing is not supplied but rather interpolated on-the-fly from 8 Hz 
satellite attitude files. 


inter-DPC exchanges and public releases, for all time-domain 
data, adding the necessary interfaces to EE files to the TOAST 
data abstraction layer. TOAST allows any analysis code to access 
any supported data set - including Planck - without worrying 
about the details of the underlying file formats or data distribu¬ 
tions on disk. Instead the path to, format of, and metadata about, 
each data file required in a particular analysis is stored in an 
XML document called a run configuration file, or “runconfig”. 
Any request by a TOAST-interfaced analysis code for data is then 
passed to TOAST to resolve and satisfy, using the information in 
the runconfig. This data abstraction also enables TOAST to cache 
heavily re-used data (such as the telescope pointing information) 
and to generate synthetic data on-the-fiy (bypassing the crippling 
10 costs associated with performing distinct time-domain simu¬ 
lation and analysis phases), and to do both of those things trans¬ 
parently to the calling analysis code. Since the general analysis 
of large, correlated data sets involves global operations on all of 
the data, requiring very massive parallelism, TOAST uses a hybrid 
MPI/OpenMP programming model to enable these operations to 
be efficiently executed on the architectures of the current largest- 
scale computing resources. All of the codes used in the various 
EEP8 pipelines are interfaced to TOAST. 

4.1. Fiducial realizations 

The EEP8 fiducial realization is generated in two steps: (1) sim¬ 
ulation of the full mission TOD for every detector; and (2) calcu¬ 
lation of maps from the various detector subsets, intervals, and 
data cuts (Eig. 11). Simulation of explicit TODs allows us to in¬ 
corporate each detector’s full beam (including its far sidelobes) 
and unique input sky (including its bandpass). As noted above, 
the fiducial realization is generated in six separate components 
- the three CMB components (lensed scalar, tensor, and non- 
Gaussian complement), two foreground realizations (with and 
without bandpass mismatch), and noise. The first five of these 
are simulated as explicit TODs and then mapped, while the noise 
is generated using the on-the-fiy approach described in the noise 
MC subsection below. 

TOD generation for any detector proceeds by first convolv¬ 
ing the appropriate sky component with the beam at every point 
in a uniformly sampled data cube of Euler angle triplets (en¬ 
coding the pointing and polarization orientation) to produce the 
“beamskyset”, and then generating the time-ordered data by in¬ 
terpolating over the beamskyset data cube to the exact point¬ 
ing and polarization orientation of each sample. Previous EEP 
simulations, including PPP6, accompanying the 2013 Planck 
data release, used the Levels software package (Reinecke et al. 
2006) to do this, using the conviqt (Prezeau & Reinecke 2010) 
and multimod tools, respectively. However, this required for¬ 
mat conversions for the input pointing data (from exchange to 
Levels) and the output time-ordered data (from Levels to ex¬ 
change), at significant 10 and disk space costs. Eor PPP8 we 
have therefore embedded the critical parts of these routines into 
a new code, lScalm2T0D, which uses TOAST to interface directly 
with exchange format data. As an additional benefit, by combin¬ 
ing the two steps into a single code the beamskyset never has 
to be written to disk, further reducing the 10 and disk space 
costs. This approach also allows us trivially to handle the LEI 
multi-component beams (subsection 3.1.3), where each beam is 
expressed at a different (^max. ^max) resolution. Previously this 
meant generating distinct TODs and adding them only after they 
had been written, tripling the 10 and disk space requirements, 
whereas now they are generated separately but combined within 
lScalm2T0D prior to being written. 
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Fig. 11. Fiducial pipeline schematic. 


All of the FFP8 fiducial maps are produced using 
Madam/TOAST, a TOAST port of the Madam generalized destrip- 
ing code (Keihanen et al. 2005, 2010), which allows for destrip- 
ing with an arbitrary baseline length, with or without a prior on 
the baseline distribution (or noise filter). Madam is used to pro¬ 
duce the official LFI maps, and its destriping parameters can be 
chosen so that it reproduces the behaviour of Polkapix, the offi¬ 
cial HFI mapmaking code. Comparison of the official maps and 
Madam/TOAST maps run using exchange data show that map- 
maker differences are negligible compared to small differences 
in pointing and (for HFI) dipole subtraction that do not impact 
the simulation. The sky components are mapped from the TODs 
generated by IScalm2T0D, while the fiducial noise is taken to be 
realization 10 000 of the noise MC (with realizations 0000-9999 
reserved for the noise MC itself). Summarizing the key differ¬ 
ences in the map making parameters for each Planck frequency: 


“ 30 GHz is destriped with 0.25 s baselines; 44 and 70 GHz are 
destriped using 1 s baselines; and 100-857 GHz are destriped 
using pointing-period baselines (30-75 min). 

- 30-70 GHz are destriped with a 1 //-shape noise prior, while 
100-857 GHz are destriped without a noise prior. 

- 30, 44, and 70 GHz have separate destriping masks, while 
100-857 GHz use the same 15 % galaxy-rpoint source mask. 

- 30-70 GHz maps are destriped using baselines derived ex¬ 
clusively from the data going into the particular map, while 
100 - 857 GHz maps are destriped using baselines derived 
from the full data set. 


4.2. Noise Monte Carlos 

The FFP8 noise MCs are generated using Madam/TOAST, ex¬ 
ploiting toast’s on-the-fiy noise simulation capability to avoid 
the 10 overhead of writing a simulated TOD to disk only to read 
it back in to map it (Fig. 14). In this implementation. Madam runs 
exactly as it would with real data, but whenever it submits a re¬ 
quest to TOAST to provide it with the an interval of the noise 
TOD, that interval is simply simulated by TOAST in accordance 
with the noise power spectral densities provided in the runconfig, 
and returned to Madam. 

For a simulation set of this size and complexity, requiring 
of the order of 10^^ random numbers over 10^^ disjoint and un¬ 
correlated intervals, care must be take with the pseudo-random 
number generation to ensure that it is fast, reliable (and specifi¬ 
cally uncorrelated), and reproducible, in particular enabling any 
process to generate any element of any subsequence on demand. 
To achieve this TOAST uses a Combined Multiple Recursive 
Generator (CMRG) (L’Ecuyer 1999) that provides more than 
sufficient period, excellent statistical robustness, and the abil¬ 
ity to skip ahead to an arbitrary point in the pseudo-random 
sequence very quickly. Rather than seed each realization sepa¬ 
rately (with the associated risk of generating overlapping sub¬ 
sequences), we then set a single seed for the entire simulation 
set, divide the resulting sequence into 2^^ disjoint subsequences, 
and assign a specific subset of subsequences to each pointing 
period of each detector for each realization, with the alloca¬ 
tion of several rather than a single subsequence allowing us 
to support multiple independent timestream components, each 
of which depends on random number generation, such as the 
HFI detectors’ uncorrelated and correlated noise components. 
Consequently when Madam requests the noise TOD for a particu¬ 
lar interval of a particular realization, the process can skip ahead 
to the appropriate location in the overall sequence; furthermore 
the random number generation itself can be threaded, with each 
thread then skipping ahead to the appropriate place within the 
subsequence. 

Madam supports MC runs by enabling the pointing and flag 
data common to all realizations to be read once, before loop¬ 
ing over the realizations that share them. However, to support 
maximum MC throughput some changes were also made to the 
standard Madam implementation: 

- once the full data for any time span are loaded, 
Madam/TOAST can independently destripe and map the full 
and two half ring data sets; 

- communication at high concurrency is more efficiently sup¬ 
ported through use of the advanced MPI collective operation 
alltoallv; 

- the map is distributed between the MPI processes based on 
pixels covered by the distributed TOD rather than using a 
simple round-robin scheme; 

- HFI full-mission destriping is supported on-the-fiy so that 
once the full mission data have been destriped, all subset 
maps can be binned as part of the same run; 

- OpenMP/MPI hybridization allows for runs at higher con¬ 
currencies without suffering from MPI communication con¬ 
gestion. 

4.3. CMB Monte Carlos 


The baseline (r = /nl = 0) fiducial channel/mission/full The FFP8 CMB MCs are generated using the FEBeCoP software 
maps for temperature and total polarization are shown in Figs. 12 package (Mitra et al. 2011), which produces beam-convolved 
and 13, respectively. maps directly in the pixel domain rather than sample-by-sample. 
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30 GHz 



44 GHz 




70 GHz 




100 GHz 



143 GHz 



217 GHz 



353 GHz 


545 GHz 


857 GHz 












liK or kJysr ^ 

Fig. 12. Left to right, the baseline FFP8 fiducial CMB, foreground and noise component temperature channel/mission/full maps 
at each frequency. Frequencies 30-353 GHz are plotted in yuK, while 545 and 857 GHz are plotted in kJy sr“^ See Fig. 16 for the 
combined temperature maps. 
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Fig. 13. Left to right, the baseline FFP8 fiducial CMB, foreground and noise component total polarization channel/mission/full maps 
at each frequency. See Fig. 17 for the combined polarization maps. 



as is done for the fiducial maps. The goal of this approach is to 
reduce the computational cost by the ratio of time-samples to 
map-pixels (i.e., the number of hits per pixel). 

As illustrated in Fig. 15, FEBeCoP proceeds in 3 steps: 

1. Given the satellite pointing and fiags and the focal plane 
(accessed through the TOAST interface), for every channel 
FEBeCoP first re-orders all of the samples in the mission 
by pixel instead of time, localizing all of the observations 
of each pixel, and writes the resulting pixel-ordered detec¬ 
tor dngles (PODA) to disk. Note that since the PODA also 
contains the detector, time-stamp, and weight of each obser¬ 


vation this is a one-time operation for each frequency, and 
does not need to be re-run for different time intervals or de¬ 
tector subsets, or for changes in the beam model or its chosen 
cut-off radius. 

2. For every time interval and detector subset to be mapped, and 
for every pixel in the map, FEBeCoP uses the PODA and the 
scanning beams to generate an effective-beam for that pixel 
which is essentially the weighted average of the discretized 
beam functions for every sample in the pixel included in the 
time interval and detector subset. The total effective-beam 
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Fig. 14. Noise Monte Carlo pipeline schematic. 



Fig. 15. CMB Monte Carlo pipeline schematic. 


array is also written to disk. Given the PODA, this is a one¬ 
time operation for any beam definition. 

3. Finally, FEBeCoP applies the effective-beam pixel-by-pixel 
to every CMB sky realization in the MC set to generate the 
corresponding beam-convolved CMB map realization. 


4.4. Validation 

Our goal for the FFP8 simulation set is that it be not only inter¬ 
nally self-consistent, but also a good representation of the real 
data. In addition to the validation steps carried out on all of the 
inputs individually and noted in their respective sections above, 
we must also validate the final outputs. A first crude level of val¬ 
idation is provided simply by visual inspection of the FFP8 and 
real Planck maps - Figs. 16 (temperature) and 17 (total polar¬ 
ization) - where the only immediately apparent difference is the 
CMB realization. 

While this is a necessary test, it is hardly sufficient, and 
the next step is to compare the angular power spectra of the 
simulated and real channel/mission/full maps. As illustrated in 
Fig. 18 and Table 5, LFI channels show excellent agreement 
across all angular scales, while HFI channels show a significant 
power deficit at almost all angular scales. Since this missing HFI 
power is not picked up in the noise estimation, it must be sky- 
synchronous (frequency bins corresponding to sky-synchronous 
signals being discarded when fitting the noise PSDs due to 
their contamination by signal residuals). This is now understood 
to be a systematic effect introduced in the HFI pre-processing 
pipeline, and we are working both to incorporate it as a system¬ 
atic component in existing simulations and to ameliorate if for 
future data releases. 

Finally, the various analyses of the FFP8 maps in conjuc- 
tion with the flight data provide powerful incidental validation. 
To date the only issues observed here are the known mismatch 
between the FFP8 and PR2-2015 cosmologies, and the miss¬ 
ing systematic component in the HFI maps. As noted above, 
the former is readily addressed by rescaling or using FFP8.1; 
however, the characterization and reproduction of the latter is an 
ongoing effort. Specific details of the consequences of this as- 
yet unresolved issue, such as its impact on null-test failures and 
p-value stability in studies of non-Gaussianity, can be found in 
the relevant papers (e.g., Planck Collaboration IX 2015; Planck 
Collaboration X 2015; Planck Collaboration XI 2015; Planck 
Collaboration XVI 2015). 


Table 5. EE and BB mission full map spectrum ratios. 


Frequency 


[GHz] EE BB 


30 . 0.998 1.000 

44 . 0.999 1.000 

70 . 0.998 1.000 

100 . 1.019 1.020 

143 . 1.033 1.029 

217. 1.017 1.007 

353 . 1.042 1.039 


This is the full mission median Q ratio (PR2-2015 / FFP8) in the 
noise-dominated range i e [2000,3000] (LFI) and i e [3000,4000] 
(HFI); as expected, FFP8.1 exhibits the same values. 


The effective-beams provide a direct connection between 
the true and observed sky, explicitly incorporating the detailed 
pointing for every detector through a linear convolution. By pro¬ 
viding the effective-beams at every pixel, FEBeCoP enables pre¬ 
cise control of systematic effects, e.g., the point-spread functions 
can be fitted at each pixel on the sky and used to determine 
point source fluxes (Planck Collaboration XXVI 2015; Planck 
Collaboration XXVII 2015). 


5. High performance computing considerations 

5.1. Simulations at NERSC 

NERSC is the US Department of Energy’s general purpose su¬ 
percomputing centre, with a new top-10 HPC system^ ^ being 

http://top56)0.org 
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Fig, 16. Pairs of channel/mission/full temperature maps, comparing FFP8 simulations (upper) and 2015 Planck data (lower). All 
maps are downgraded to Aside = 256. Visible differences away from the galactic plane are between the actual and simulated CMB 
sky. 
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Fig. 17. Pairs of channel/mission/full total polarization maps, comparing FFP8 simulations (upper) and 2015 Planck (lower) data. 
All maps are downgraded to Aside = 256. 
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Fig. 18. Comparison of the BB spectra of the channel/mission/full simulated and flight data maps using pseudo-spectra computed on 
15% of the sky with the Galaxy and point sources masked, showing excellent agreement at 70 GHz but a few percent discrepancy at 
almost all angular scales at 100 GHz. Left: 70 GHz full map 55-spectra. Right: 100 GHz full map 55-spectra. Top: Linear horizontal 
axis shows small scale behaviour. Bottom: Logarithmic horizontal axis shows large scale behaviour. 


installed about every 3 years. For the last 15 years it has pro¬ 
vided annual allocations of cycles and storage to the Planck col¬ 
laboration, with a unique agreement between NASA and DOE 
guaranteeing Planck access to NERSC through the mission life¬ 
time (Borrill et al. 2014). The US Planck collaboration has also 
purchased exceptional levels of service at NERSC, including 
240 TB of globally accesssible disk space and one cabinet (640 
cores) of the 9 984-core IBM iDataplex Carver. In 2014, Planck 
was awarded about 75 000000 NERSC CPU-hours, which were 
shared by the 170 Planck data analysts who have requested 
NERSC accounts. 

Eor EEP8, the flducial realization, all HEI noise MCs, 
the first 1 000 LEI channel noise MCs, and all CMB MCs 
were performed on NERSC’s 133 824-core Cray XC30 Edison. 
Subsequently, all of the EEP8.1 CMB MCs were generated there 
too. 

5.1.1. The fiducial realization 

The EEP8 flducial foreground and CMB components were gen¬ 
erated via explicit timelines and the noise via an implicit time¬ 


line, as detailed in the noise MC section below. Eor the ex¬ 
plicit cases each frequency used a separate job script, allow¬ 
ing multiple frequencies to be executed simultaneously without 
concerns about load-balancing or synchronization. Within each 
frequency-specific job each component was run in turn, allow¬ 
ing the same target TOD files to be re-used (overwritten) to save 
disk space. Within each component, each detector was run in 
turn to prevent overloading the file system (since the exchange 
format TODs are stored by frequency). Once a full set of TODs 
had been generated for a given component at a given frequency, 
the full set of 1134 maps was produced before moving on to the 
next component. The concurrencies and runtimes for each fre¬ 
quency are listed in Table 6. With five components generated 
via explicit timelines and one via on-the-fly simulation within 
the mapmaking, the total cost was approximately 250 000 CPU- 
hours. This same set-up also allowed us to generate each individ¬ 
ual foreground component to cross-check the component sepa¬ 
ration pipelines. 


22 











Planck Collaboration: Full Focal Plane Simulations 


Table 6. Computational cost of generating one component of the 
FFP8 fiducial realization for each frequency. TOD generation 
includes all of the detectors, and mapmaking all of the maps, at 
that frequency. Note the much higher cost of LFI mapmaking, 
due to its shorter baselines and run-baseline destriping. 


Frequency Concurrency TODs Maps Cost 

[GHz] (Cores) [min] [min] [CPU-hrs] 


30 . 3 600 12 20 1920 

44 . 3 600 36 15 3 060 

70 . 3 600 108 48 9 360 

100 . 4 800 64 3 5 360 

143 . 4 800 88 5 7440 

217 . 4 800 96 5 8 080 

353 . 4 800 96 5 8 080 

545 . 4 800 18 2 1600 

857 . 4 800 24 2 2080 


5.1.2. The noise Monte Carlos 

For the noise MCs, we noted that mapmaking runs most ef¬ 
ficiently at the lowest concurrency with sufficient memory to 
hold all of the data at a given frequency. Additional parallelism 
was then achieved by running multiple instances of a single fre¬ 
quency simultaneously. In practice, for each instrument and fre¬ 
quency we first ran a job containing 10 instances, each gener¬ 
ating 10^ realizations of all map-types, to produce the first 10^ 
realizations, and then, for HFI only (for LFI, see Sect. 5.2), ran 
a second job running nine instances each generating 10^ real¬ 
izations of only the channel/mission/full maps to fill out the full 
10"^ realizations. For LFI, where every data subset is individually 
destriped, the instances were further subdivided to keep the data 
per core constant, so that the same number of cores would pro¬ 
cess one channel/mission, two channel/half-mission, four chan¬ 
nel/year, or eight channel/survey maps, or (for the 70 GHz data) 
three detset/mission, six detset/half-mission, twelve detset/year 
or twenty-four detset/survey maps. Table 7 shows the cost per 
realization of all the maps at one frequency (for realizations 0- 
999) and of the HFI channel/mission/full maps alone (for real¬ 
izations 1000 - 9999). The overall cost of these runs was ap¬ 
proximately 3 000000 CPU-hours for the 10^ realizations of all 
maps and 5M CPU-hours for the 9 x 10^ additional HFI chan¬ 
nel/mission/full maps. 

Table 7. Computational cost of generating one realization of the 
FFP8 noise MC for each frequency, both for the full set of maps 
and for the HFI channel/mission/full (CMF) maps alone. 


All maps CMF maps 

Frequency Concurrency _ _ 

[GHz] (Cores) [min] [CPU-hrs] [min] [CPU-hrs] 


30 . 2 304 10 384 . 

44 . 3 456 8 460 . 

70 . 4 608 13 949 . 

100 . 2 640 3 176 2 88 

143 . 3 600 3 180 2 120 

217 . 3 960 3 198 2 132 

353 . 3 960 3 198 2 132 

545 . 600 4 40 2 20 

857 . 600 6 60 3 30 


5.1.3. The CMB Monte Carlos 

For the CMB MCs, the computational cost depends on the num¬ 
bers of time-samples TVt, map-pixels TVp, beam-pixels A/b, and 
MC realizations Generating the PODA is essentially a re¬ 
ordering of the time-domain data, scaling as Nu generating the 
effective-beam is a calculation of beam-pixel weights for ev¬ 
ery observation, scaling as A(tA/b; and generating the beam- 
convolved maps is an effective-beam-matrix/map-vector multi¬ 
plication for each realization, scaling as A/^mcA/pWb. Table 8 
shows the values of these parameters for each Planck channel, 
indicating that for FFP8 the effective-beam convolution over the 
CMB skies has, for the first time, become the dominant cost. 
By comparison, for FFP6, although the numbers of samples 
were a factor of 2 (HFI) to 4 (LFI) smaller, the smaller MC 
set (with only 10^ realizations) and smaller numbers of pixels 
(temperature only) made the effective-beam calculation domi¬ 
nate. Running the 46 distinct lO'^-realization CMB MC map sets 
included in FFP8 then required some 8M CPU-hours in total. 


Table 8. The numbers of samples (TVt), map-pixels (A/p), 
beam-pixels (A/b) and the derived computational costs of the 
FEBeCoP effective-beam calculation (TVt A/b) and MC convolu¬ 
tion (A/^mc >/Vp M)) ior each channel in FFP8. The number of 
effective-beam pixels at each frequency depends on the HEALPix 
resolution of the map (A/side = 2048 and 1024 for HFI and LFI, 
respectively) and the cut-off radius of the scanning beam (113.6', 
79', 52' and 100' for 30 GHz, 44 GHz, 70 GHz and HFI, respec¬ 
tively). 


Frequency 

K 



KK 

TVmc TVp TVb 

30 ... . 

1.5x10*0 

3.8x100 

3.5x10^ 

5.2x10*’ 

1.3x10*’ 

44 ... . 

3.2x10*0 

3.8x10** 

1.7x10’ 

5.4x10*’ 

6.4x10*“* 

70 ... . 

1 .1x10** 

3.8x10** 

7.3x10’ 

7.9x10*’ 

2 .8x10*“* 

100 ... . 

7.5x10*0 

1.5x10’ 

l.lxlO-* 

8.3x10*“* 

1.7x10*0 

143 ... . 

1 .0x10** 

1.5x10’ 

1 .1x10“* 

1 .1x10*’ 

1.7x10*0 

217 ... . 

1 .1x10** 

1.5x10’ 

1 .1x10“* 

1 .2x10*’ 

1.7x10*0 

353 ... . 

1 .1x10** 

1.5x10’ 

1 .1x10“* 

1 .2x10*’ 

1.7x10*0 

545 ... . 

3.2x10*0 

5.0x100 

1 .1x10“* 

3.6x10*“* 

5.5x10*’ 

857 ... . 

4.3x10*0 

5.0xl00 

1 .1x10“* 

4.8x10*“* 

5.5x10*’ 


Originally FEBeCoP was designed to optimize MC simu¬ 
lations with well-localized beams (and hence small Wb), but 
it is clear that with FFP8 we have moved out of that regime 
and it is time to re-evaluate the comparison of the computa¬ 
tional cost and accuracy of the effective-beam and explicit time- 
domain approaches. Since these scale as TVp TVb and Wt, respec¬ 
tively, the cross-over point is when the number of pixels in the 
effective-beam is comparable to the number of hits per map- 
pixel. Possible ways to reduce TVb for FEBeCoP would be to use 
multi-component beams with each component pixelized at a dif¬ 
ferent resolution, or to restrict the main beam to the pixel domain 
and process the intermediate and far sidelobes in the spherical 
harmonic domain. Looking forward, this will be an active area of 
research and development. We note, though, that next-generation 
5-mode experiments are specifically designed to have very many 
more hits per pixel than Planck in order to achieve the necessary 
signal-to-noise ratio to meet their science goals, shifting the bal¬ 
ance back towards pixel-domain calculations. 
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5.2. Simulations at CSC 

A significant fraction of the LFI noise MCs were performed 
at the CSC-IT Center for Science in Finland, using their Cray 
XC30 supercomputer Sisu with 40512 computing cores. These 
runs were performed using the same software stack as was used 
at NERSC, with cross-checks confirming that the two systems 
produce results identical to numerical precision. 

The CSC runs included the channel/mission/full and half¬ 
ring maps for realizations 1000-9 999, totalling 81000 maps 
and using a total of 4 million CPU-hours. These were run in 
about a dozen different jobs, gradually increasing the number of 
instances of the code running simultaneously and the number of 
cores used. The largest jobs for each frequency were as follows: 

- for 30 GHz, 10 simultaneous instances generating 5 000 re¬ 
alizations in all, running on 23 040 cores for 17.8 hours for a 
total of 410 000 CPU-hours; 

“ for 44 GHz, 10 simultaneous instances generating 5 000 re¬ 
alizations in all, running on 11 520 cores for 21.8 hours for a 
total of 251 000 CPU-hours; 

- for 70 GHz, 10 simultaneous instances generating 2757 re¬ 
alizations in all, running on 36 000 cores for 24 hours^^ for a 
total of 864 000 core-hours. 

In addition, we produced a smaller number of realizations 
for the 70 GHz detector subsets, including 1 000 realizations of 
the mission and 300 realizations each of the two half-missions, 
four years, and eight surveys. With three detector subsets and 
three maps (full and half-ring) for every realization, this resulted 
in 46 800 maps and took a total of 200 000 CPU-hours. 

In total, 127 800 maps were produced at CSC, using 
4 200000 CPU-hours and requiring a total of 17.7 TB of disk 
space; these were then transferred to NERSC for wider distribu¬ 
tion. 


6. Results 

The core EEP 8 simulation set consists of: 

- 90720 fiducial maps, comprising 1 134 maps of each of six 
components (three CMB, two foreground, one noise), then 
combined to produce two sets of total maps (with and with¬ 
out bandpass mismatch), all for five different values of the 
(^5 /nl) parameter pair. 

“ 671 400 noise MC maps, comprising 10"^ realizations of all 
channel/mission/full maps, at least 10 ^ realizations of all 
other full maps, and at least 10 ^ realizations of all half-ring 
maps. 

- 460000 CMB MC maps, comprising 10"^ chan¬ 
nel/mission/full and channel/half-mis sion/full maps and 
selected detset/mission/full maps. 

while the EEP8.1 simulation set currently consists of: 

- 18 144 fiducial maps, comprising 1 134 maps of each of three 
CMB components, then combined with the EEP 8 foreground 
and noise maps to produce two sets of total maps (with and 
without bandpass mismatch), restricted the the baseline (r = 
/nl = 0 ) case. 

- 27 000 CMB MC maps, comprising 1 000 chan¬ 

nel/mission/full and channel/half-mission/full maps. 

The maximum runtime allowed at CSC for jobs of this size. 


6.1. Uses of FFP8 

As well as supporting the validation and verification of all of 
the Planck analysis codes, EEP8 is widely used in the analysis 
of Planck data. 

Fiducial realizations are used to: 

- test how various choices of the smoothing kernel af¬ 
fect signal aliasing in low resolution maps (LEI) (Planck 
Collaboration VI 2015); 

- demonstrate inaccuracies in the measurements of the instru¬ 
mental bandpasses (Planck Collaboration X 2015); 

- measure and subtract the LEI sidelobe pick-up (Planck 
Collaboration II 2015); 

- estimate calibration and bandpass measurement uncertain¬ 
ties in component separation (Planck Collaboration X 2015); 

- de-bias and measure uncertainty in lensing reconstruction 
(Planck Collaboration XV 2015); 

- assess SZ catalog reliability (Planck Collaboration XXVII 
2015); 

- estimate noise for mapping the thermal Sunyaev-Zeldovich 
(tSZ) effect (Planck Collaboration XXII 2015). 

Noise MCs are used to: 

- measure the accuracy of the LEI pixel-pixel noise covariance 
matrix construction (Planck Collaboration VI 2015); and 

- inform low-^ data selection by measuring the distribution 
of low-multipole amplitudes due to noise as well as con¬ 
firming the anomalous nature of Surveys 2 and 4 (Planck 
Collaboration II 2015); 

CMB MCs are used to: 

- produce beam transfer functions for each map type listed in 
Table 1 by averaging over the corresponding MC realiza¬ 
tions (Planck Collaboration IV 2015; Planck Collaboration 
VII 2015). 

In addition, the noise and CMB MC maps are com¬ 
bined using the same component separation pipelines as the 
real and fiducial data. The resulting MC sets are then used 
to support the various CMB isotropy and statistics analyses 
(Planck Collaboration XVI 2015) and non-Gaussianity estima¬ 
tors (Planck Collaboration XVII 2015). 

Various ancilliary data sets were also produced to address the 
specific requirements of particular analyses. The most notable of 
these were: 

- nine channel/mission/full maps for each of the ten individual 
foreground components, to validate the component separa¬ 
tion analyses; 

- 10^ channel/mission/full of the scalar CMB MCs with no 
Doppler boosting applied, to validate the BipoSH reconstruc¬ 
tion of Doppler boosting (Planck Collaboration XVI 2015). 

The capabilities developed for PPP8 also allowed for many 
other tests of data and analysis systematics, including running 
the mapmaking: 

- with both mission- and run-baseline destriping to measure 
correlations induced by mission-baseline destriping in maps 
spanning different time intervals or including disjoint detec¬ 
tor subsets, leading to the acceptance of this approach by 
HEI (where these are generally small) and its rejection by 
LEI (where they are significant); 
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- with both temperature-only and temperature-rpolarization 
destriping to demonstrate the loss of polarization power in¬ 
duced by the former, leading to the adoption of the latter by 
both HFI and LFI; 

- using pointing that includes the expected uncertainty in its 
reconstruction to test the impact of pointing reconstruction 
error for LFI maps; 

“ using polarization angles that include the expected uncer¬ 
tainty in their determination to test the impact of polarization 
angle error for LFI maps. 

6.2. Known issues with FFP8 

Generating a simulation set of the size and complexity of FFP8 
requires the coordination of a large number of people and pro¬ 
cesses. In addition to the cosmology mismatch and HFI system¬ 
atic residual issues noted above, there are a few errors that are 
known to have been inadvertently introduced into FFP8 during 
its execution. The most significant of these derives from an ap¬ 
proximation in the way the PSM applies the detector bandpasses 
to the input skies, which results in effective bandpasses that dif¬ 
fer from the inputs by up to 2% in the submillimetre chan¬ 
nels. This effect has been incorporated into the fiducial maps, 
but when the submillimetre CMB MC maps are used with the 
fiducial realization they must be rescaled by factors of 0.9976 at 
545 GHz and 0.9828 at 857 GHz. 

The other known issues are not expected to have any signifi¬ 
cant impact, but we list them here for completeness. 

- During the export of the TOD from the HFI DPC the sample 
time stamps were inadvertently truncated, resulting in point¬ 
ing errors of up to 0.'3. 

- The channel-average bandpasses generated by the PSM used 
preliminary detector weights, resulting in a tiny discrepancy 
between the bandpass-mismatch and no-bandpass-mismatch 
cases. 

- The FFP8 mapmaking used a slightly more rigorous pixel- 
rejection criterion than the HFI DPC, resulting in some FFP8 
maps having somewhat more missing pixels than the real 
data. Since these are the noisiest pixels, their naive inclu¬ 
sion in any map-based analysis can have a disproportionate 
effect at high i. 

6.3. Community access to FFP8 

As was the case with FFP6 for the first Planck data release, 
FFP8/FFP8.1 is being made available to the scientific commu¬ 
nity. A simulation set of this size requires significant computa¬ 
tional resources to store and analyse, so once again we are mak¬ 
ing it available at NERSC, where there is also an allocation of 
compute cycles available to support its exploitation. Full details 
on this, including how to obtain a NERSC account, can be found 
at http://crd.Ibl.gov/cmb-data. 

7. Conclusions 

Simulations play a critical role in the analysis of CMB data 
sets, most particularly in the validation and verification of anal¬ 
ysis pipelines and the debiasing and uncertainty quantification 
of analysis results. The former requires a simulation that is as 
true to the mission as possible, while the latter requires enor¬ 
mous Monte Carlo sets of mission realizations. FFP8 provides 
the 2015 Planck data release with both. The fiducial realization 
takes the most realistic model of the microwave sky available to 


us and simulates the full mission observation of that sky for ev¬ 
ery Planck detector, incorporating our best estimates of that de¬ 
tector’s beam, bandpass, and time-varying noise properties. The 
noise and CMB MC sets include 10"^ realizations of the most 
important maps (namely those that incorporate all of the data 
at a given frequency), together with 10^-10"^ realizations of all 
of the other maps that include subsets of the observations and 
detectors at each frequency. The FFP8.1 CMB maps provide a 
smaller sample of fiducial and MC CMB skies consistent with 
the PR2-2015 cosmology. With around 1.25 million maps in to¬ 
tal, and requiring some 25 million CPU-hours at multiple HPC 
centres, this is by far the largest CMB simulation set ever pro¬ 
duced. 
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Appendix A: Measuring band-pass mismatch 
residuais through spurious-component maps 

We measure bandpass mismatch by projecting the mismatch into 
extra degrees of freedom called “spurious” maps. The mismatch 
has a specific signature: it is a sky-synchronous difference be¬ 
tween the detectors at a given frequency that does not modu¬ 
late under rotation (unlike polarization). We have extended the 
method to measure the mismatch between individual detectors 
and the noise-weighted frequency average. Conventionally, the 
mismatch is only resolved between two polarization-orthogonal 
detectors that have common pointing (i.e., are co-located on the 
focal plane). 

Mathematically we can formulate this as an extension of the 
pointing model. We start from the conventional polarized detec¬ 
tor response: 

= Ip + rf* (Qp cos 2{at + if/^) + Up sin 2{at + j/^)) + fit, 

(A.l) 

where detector D TOD sample, at time f is a linear combi¬ 
nation of the observed intensity sky, /, at pixel p and linearly 
polarized (2, U) sky that is modulated by the polarization effi¬ 
ciency, T]^, the polarization sensitive angle of the detector, 
and the orientation of the detector at. In this model, the sky sig¬ 
nal includes the azimuthally symmetric part of the beam, and all 
instrumental noise and systematics are contained in the noise, n. 
The full data vector can be written concisely as 

d = Pm -f n, (A.2) 


where P is the pointing matrix that scans the smoothed sky map, 
m, into the time domain according to the instrument scanning 
pattern. 

We extend the model by assuming that the detectors are ac¬ 
tually observing sky emission that is convolved with the individ¬ 
ual detector bandpass, causing each detector to view a slightly 
different sky, The difference between an effective frequency 
average, /, and the sky seen by a particular detector, /^, is the 
spurious (or bandpass mismatch) map, = I - . We can 

solve for the frequency-averaged (/, 2, U) sky by introducing 
a hierarchy of spurious maps that captures the mismatch first 
within pairs of detectors (horns), then between pairs of horns 
(quadruplets), and ultimately between larger aggregates. We call 
all these combinations “virtual horns”. We show in Table A.l 
how the spurious maps can be assigned for the Planck 70 GHz 
channel. The specific hierarchy of virtual horns is not impor¬ 
tant: it is only required that for each pair of detectors there exists 
exactly one spurious map (virtual or physical horn) for which 
they have opposite pointing weights. The total bandpass mis¬ 
match map between any detector and the frequency average can 
be produced as a linear combination of the spurious maps with 
the appropriate pointing weights. 

The spurious maps extend the pointing model in Eq. A. 1 to 

+ (A.3) 

i 

where the sum is over all spurious maps that correspond to phys¬ 
ical and virtual horns listed in Table A.l and the associated 
weights are -hi, -1, or 0. These extra pointing weights are in¬ 
cluded in additional columns in the pointing matrix, P, and a 
full map with all the spurious components can be binned using 
the standard mapmaking equation: 

m = (P'^N-' P)”‘ P^N-irf, (A.4) 

where N is some model of the detector noise covariance, fre¬ 
quently a diagonal matrix with the detector white noise vari¬ 
ances. In practice, the additional degrees of freedom in the spu¬ 
rious maps can only be resolved for pixels that have sufficient 
coverage in crossing angles to break the degeneracy between po¬ 
larization and bandpass mismatch. Even then, the resulting noise 
levels in the /, Q, and U maps are much higher than in the stan¬ 
dard mapmaking case, making the IQUS solutions only useful 
for low-resolution studies and checking other bandpass leakage 
correction schemes that rely on prior information such as detec¬ 
tor bandpasses and models of the sky emission. 
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Table A.l. Spurious map pointing weights 


Detector 
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Homs 1-6 are the physical horns that comprise two polarization-orthogonal arms. 

Homs 7-9 are the “virtual horns”; each consists of two physical horns, LFI18-hLFI19, LFI20-hLFI21, and LFI22-hLFI23. 
Homs 10-11 expand the hierarchy to include mismatch between the virtual horns. 
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